Statistical laws of stick-slip friction at mesoscale

Friction between two rough solid surfaces often involves local stick-slip events occurring at different locations of the contact interface. If the apparent contact area is large, multiple local slips may take place simultaneously and the total frictional force is a sum of the pinning forces imposed by many asperities on the interface. Here, we report a systematic study of stick-slip friction over a mesoscale contact area using a hanging-beam lateral atomic-force-microscope, which is capable of resolving frictional force fluctuations generated by individual slip events and measuring their statistical properties at the single-slip resolution. The measured probability density functions (PDFs) of the slip length δxs, the maximal force Fc needed to trigger the local slips, and the local force gradient \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}^{{\prime} }$$\end{document}k′ of the asperity-induced pinning force field provide a comprehensive statistical description of stick-slip friction that is often associated with the avalanche dynamics at a critical state. In particular, the measured PDF of δxs obeys a power law distribution and the power-law exponent is explained by a new theoretical model for the under-damped spring-block motion under a Brownian-correlated pinning force field. This model provides a long-sought physical mechanism for the avalanche dynamics in stick-slip friction at mesoscale.


A. Characterization of rough and smooth substrates
An ultra-fine silicon carbide sandpaper taped on a flat glass slide is used as the rough substrate for the friction measurement.The sandpaper has an average grain size of 100 nm, as shown by the AFM topographical image in Fig. 1(c) of the main text.Here we provide additional characterizations of the sandpaper surface.Supplementary Fig. 1(a) shows an AFM topographical image of the sandpaper over an area of 50 × 50 µm 2 , from which we find the RMS roughness of the sandpaper surface to be 0.17 µm (over an area of 50×50 µm 2 ).It is also seen that the roughness of the sandpaper is quite uniform with very few large-sized pits (> 10 µm).This uniform roughness allows the scanning probe to move steadily without being fully trapped during the scan.Supplementary Fig. 1(b) shows the roughness height profile h(x) along the red Supplementary Fig. 2. Characterization of surface roughness of the sandpaper.(a) Measured PDF P (h) of surface height h.The mean value of h is set to zero.The red line shows a Gaussian fit to the data points with a standard deviation of 0.17 µm.(b) Measured power spectral density (PSD) functions S(q) of the roughness height profile h(x) as a function of the wave-vector q = 2π/δ (bottom axis).The top axis indicates the corresponding wave-length δ.The red circles are obtained from the AFM topographical image shown in Fig. 1(d) of the main text (10 × 10 µm 2 ) with a spatial resolution of 0.039 µm.The black squares are obtained from the AFM topographical image shown in Fig. 1(a) (50 × 50 µm 2 ) with a spatial resolution of 0.049 µm.The dashed lines indicate the two power-law regions with the numbers indicating the power law exponents.line marked in Fig. 1(a).The abrupt variations in the obtained h(x) indicate sharp protrusion of the grain clusters on the sandpaper.These grain clusters (or asperities) will exert strong friction (or pinning) to the scanning probe.
We now discuss the statistical properties of the roughness height profile h(x) of the sandpaper.Supplementary Fig. 2(a) shows the probability distribution function (PDF) P (h) of the surface height h, where the mean value of h is set to zero and h is sampled over an area of 50 × 50 µm 2 .It is found that the measured P (h) is well described by a Gaussian distribution with a standard deviation 0.17 µm, which is the RMS roughness.Supplementary Fig. 2(b) shows the measured power spectral density (PSD) functions S(q) of the roughness height profile h(x).Because the sandpaper surface is isotropic, herein we only consider the one-dimensional PSD.It is seen that the two PSD functions obtained at two different spatial resolutions overlap well in the common range of the wave-vector q.The functional form of the measured S(q) over the entire q-range (10 5 m −1 q 10 8 m −1 ) can be described by two regimes of approximately power-law scaling, S(q) ∼ q −α , together with a narrow intermediate cross-over region centered around q ∼ 10 7 m −1 .The power-law exponent α in the small-q regime (or long-wavelength regime with δ 2 µm) is α ≃ 1.6 (left dashed line) whereas the value of α in the large-q regime (or short-wavelength regime with δ 0.3 µm) is α ≃ 3.5 (right dashed line).While the cross-over region has a finite width, we find that the two power-law fits to the data meet at a cut-off length ℓ c ≃ 0.34 µm.Similar cross-over behaviors from one power-law-like regime to another were also observed for other rough surfaces [1][2][3].
A power-law-like regime in the measured PSD suggests that the rough surface has some self-similar features and such a surface is referred to as a fractal.By examining the AFM images (Fig. 1(a) and Fig. 1(d) of the main text) of the sandpaper surface, we find that the powerlaw-like regime at small δ (from ∼ 0.1 µm to ∼ 0.4 µm) can be attributed to the individual gains (with an average grain size 0.1 µm) and their small clusters (such as dimers, trimers and tetramers).A similar conclusion was also drawn in a previous study [3].On the other hand, the power-law-like regime at large δ (from ∼2 µm to ∼30 µm) is caused primarily by the random distribution of larger clusters (aggregation of small clusters), which are associated with the bright regions (or bumps) and dark regions (or depressions) in the AFM topographical images.The cut-off length ℓ c that separates the two powerlaw-like regimes, therefore, can be viewed as the lower cutoff of the large cluster (or domain) sizes.Supplementary Fig. 2(b) thus demonstrates that the sandpaper has various asperities (or micro-contacts) with a broad range of sizes.
For our system, the maximum micro-contact size is limited by the dimension of the apparent contact area along the scanning direction, which is ∼ 4 µm for the quasi-1D probe and ∼ 12 µm for the 2D probe.The minimum micro-contact size, on the other hand, is the size of the grain particles on the sandpaper.Consequently, when the AFM probe scans over the sandpaper, it will experience various asperities (or microscopic contacts) with a broad range of scales from the grain size of ∼ 0.1 µm to the probe size up to ∼ 12 µm.Within the two length limits, the measured PSD (as a function of wavelength δ) has two (short-ranged) power-law-like regimes, which nevertheless cannot be described by a single fractal.
As mentioned in the main text, the main purpose of using two contact materials with a large difference in hardness is to increase the compliance of the contact so that the number of asperities (or micro-contacts) with different random configurations (e.g., different roughness heights and sizes) is maximized at the contact without introducing large plastic deformations and wear.More specifically, the rigid sandpaper surface used here serves as a base to provide a large number of asperities with a broad range of randomness, whereas the probe surface is made of a softer material and is flat so that it is able to have a close contact with the rough surface of the sandpaper at all length scales and feel the underlying pinning force field, when a small normal load is applied.If the probe surface were made of another rigid sandpaper, for example, the contact between the two hard surfaces is dominated primarily by a small number of the most prominent asperities on each surface.In this case, it would be difficult to have an effective local contact under a regular normal load that can feel many asperities with a broad range of randomness.
To measure the static spring constant k 0 of the scanning probe in contact with a non-slip substrate, we use a smooth silicon wafer surface coated with a thin layer of gold as the non-slip substrate.Supplementary Fig. 3 shows an AFM topographical image of the gold-coated silicon wafer surface, which has a RMS roughness of 1.0 nm over an area of 10 µm×10 µm.

B. Calibration of lateral force
When a lateral (or normal) force is applied to the AFM tip, the cantilever will bend laterally (or vertically) causing a displacement of the laser spot on the photodetector of the AFM.This displacement is recorded by the AFM as a voltage signal proportional to the magnitude of the applied force.To convert this voltage signal to the actual force, one needs to calibrate two proportionality constants.One is the optical sensitivity of the detector, which is the displacement of the tip per voltage change; the other is the spring constant of the probe, which is the force per tip displacement.For the vertical bending mode, the optical sensitivity is extracted from the slope of the voltage-displacement curve measured when the AFM cantilever is pressed against a flat hard surface.The vertical spring constant is obtained using the thermal power spectral density (PSD) method [4,5].Using this method, we find the normal spring constant of the AFM cantilevers used in this experiment to be k n ≃ 0.7 N/m.
For the lateral bending mode used in the friction measurement, however, the optical sensitivity usually cannot be obtained using the same method as that for the vertical bending mode because of the undetermined stiffness at the solid interface [6][7][8][9].Furthermore, the thermal PSD method is not applicable to the lateral bending mode.As a result, the lateral spring constant of the cantilever together with the tip is usually estimated based on the sample dimensions.
In this experiment, we make a direct calibration of the lateral force sensitivity, the coefficient that relates the measured voltage signal to the actual lateral force, using a glass fiber as a force sensor.The glass fiber is pulled from a glass rod of 1 mm in diameter using a pipette puller (P-97, Sutter Instrument) and has the dimensions of ∼ 10 µm in diameter and ∼ 5 mm in length.The fiber is linked to a segment of unstretched glass rod (1 mm in diameter) as a handle for convenience of manipulation.The spring constant of the glass fiber is calibrated by vertically pressing an AFM cantilever of known spring constant (10.9 N/m) against the free end of the horizontally aligned glass fiber, as shown in Supplementary Fig. 4(a).The slope of the measured force-displacement curve is the effective spring constant of the joint cantilever-fiber system, which is composed of two springs connected in series.With the known spring constant of the cantilever, we obtain the spring constant of the glass fiber to be ≃ 0.029 N/m, which is very small compared with that of the cantilever.
A similar method is used to calibrate the lateral force sensitivity of the scanning probe.As shown in Supplementary Fig. 4(b), the scanning probe is mounted on the AFM cantilever holder and the glass fiber is mounted horizontally on the AFM scanner with its long axis in parallel with the cantilever arm of the scanning probe.Using the controllers of the AFM cantilever holder and scanner, the glass fiber is moved gradually toward the scanning probe until its tip touches the scanning probe.We then obtain the lateral force-displacement curve, where the measured force is in unit of volt and the displacement is the total deflection of the glass fiber and the scanning probe.Because the spring constant of the glass fiber is much smaller than the lateral spring constant of the scanning probe (∼ 1 N/m), the deflection of the scanning probe is negligible.Therefore, by multiplying the displacement by the spring constant of the glass fiber, we convert the force-displacement curve to force (in unit of volt)-force (in unit of Newton) curve and the slope is the proportionality constant between the measured voltage signal and the actual lateral force.

C. Effective lateral spring constant of the scanning probe
The effective lateral spring constant k e of the scanning probe as a whole is determined by the contributions from its three components, namely, the horizontal cantilever (k c ), the vertical hanging beam (k b ), and the beam tip with UV cured glue in contact with the substrate (k t ).Because the three springs are connected in series, we have [6] 1 The spring constant of the individual components can be determined by their dimensions and material properties.The lateral spring constant of a rectangular AFM cantilever with a hanging tip is given by [6,7] where G is the shear modulus of the cantilever (G ≃ 51 GPa for silicon [10]), and the cantilever has the dimensions of length L, width W , thickness T and tip length ℓ (or the hanging beam length ℓ).For the probe used in the experiment, we have L ≃ 300 µm, W ≃ 35 µm, T ≃ 2 µm, ℓ ≃ 90 µm, and the resulting k c ≃ 1.9 N/m.The lateral spring constant of the hanging beam is actually the normal spring constant when used as an AFM cantilever and is given by [6,7] where E is the Young's modulus of the hanging beam (E ≃ 169 GPa for silicon [10]), and the beam has the dimensions of length ℓ, width w, and thickness t.With ℓ ≃ 90 µm, w ≃ 35 µm, and t ≃ 2 µm, we have k b ≃ 16 N/m.The lateral spring constant of the hanging beam tip in contact with the substrate can be estimated by assuming that the UV cured glue layer and the embedded hanging beam tip functions as a rubber block with an effective shear modulus G t .In this case, we have where A r is the real contact area and ℓ is the hanging beam length.When the contact between the hanging beam tip and the substrate is "perfect" such that the real contact area A r is close to the apparent contact area A a , k t reaches its maximal value k max t .Because the shear modulus of the substrate (either the sandpaper or the Au-coated silicon wafer) is much larger than that of the nanoparticle-embedded UV cured glue, the effective shear modulus G t should be similar to the shear modulus of the UV cured glue, which is estimated as  1).The estimated k max e is in the same order as the measured static spring constants, k 0 ≃ 0.66 N/m for the quasi-1D probe and k 0 ≃ 0.48 N/m for the 2D probe.The difference between the two values may result from an overestimation of the real contact area A max r in Eq. (4) due to the surface roughness of the hanging beam tip with a nanoparticle-glue coating.

D. AFM operation
In the experiment, an AFM (MFP-3D, Asylum Research Inc.) is used for the friction measurement.The AFM operates in the lateral force microscopy mode, in which the scanner drives the substrate (the sandpaper) moving in the direction perpendicular to the cantilever.Before the scan starts, the scanning probe is pressed against the sandpaper under a normal load N .The probe scans on the sandpaper at a constant speed U in one direction for a distance of 100 µm followed by a re-scan along the same track in the reversed direction.During the scanning, the displacement and lateral force signals are recorded at a sampling rate 1 kHz.
In the meso-scale friction experiment, one may examine the level of interfacial wear by checking the reproducibility of the frictional force trajectories F (x) obtained when the scanning probe repeats its scans over the same region of the sandpaper surface and under the same condition.Supplementary Fig. 5 shows three such repeated force trajectories.The three force trajectories overlap with each other in most positions, suggesting that x (mm) Supplementary Fig. 5. Reproducibility of the frictional force measurements.The black, red and blue curves are obtained sequentially along the same track on the sandpaper under the same scanning conditions: normal load N = 400 nN and scanning speed U = 100 nm/s.Three scans are made and the total travel distance for each scan is 200 µm.
the contact geometry of the interface remains essentially unchanged during the repeated scans and the interfacial wear is not significant.There are a few locations where the three curves show some deviations.These deviations occur at the locations where a large slip is replaced by two or more smaller ones (such as those shown at x = 54.2 µm and x = 55.3 µm) or small-sized slips take place repeatedly (such as those shown at x = 54.8 µm and x = 57.6 µm).The deviations may be caused by the stochastic slip dynamics or by some minor plastic deformation of small asperities.With these results, we conclude that by ignoring the plastic deformation of small asperities (if any), the interfacial wear in our system is negligible.To obtain the probability density function with enough independent slip events, we repeat the force measurements five times along different tracks on the sandpaper under the same scanning condition.A total of ∼2500 slip events are collected over a total scanning distance of 500 µm under each scanning condition.
A common issue associated with the frictional force measurement using AFM is that the normal load N and the lateral force F are partially coupled because the torsional bending of the cantilever involves both the vertical and lateral components.Therefore, fluctuations of the lateral force F will result in fluctuations of the normal load N .To keep N constant, the deflection feedback control (for normal load) is usually used, which is very helpful for measurement of the frictional force over a flat surface where fluctuations of the frictional force are small.However, the feedback will introduce artificial discontinuities to the measured frictional force over rough surfaces.Suppression of the deflection feedback is, therefore, needed for the AFM study of stick-slip friction.In our experiment, we set the gain of the deflection feedback to a small value (∼0.001) so that the feedback is close to an off state.This is equivalent to set the lateral scan at a constant "indentation".This setting has an advantage over that of completely turning off the feedback control in that the normal load can still be tuned slowly in case there is an incline on the substrate.Supplementary Fig. 6 shows how the deflection feedback gain affects the measured frictional force F .When the gain is set at 0.001 (red curve), the measured frictional force reveals a clear zigzag pattern for the stick-slip motion of the hanging bean probe.When the gain is set at 0.1 (blue curve), however, the measured F shows many irregular discontinuities, as indicated in Supplementary Fig. 6.These irregular spikes occur when the probe is at the stick state as a result of the fast change of the normal load.When the scanning probe slides against the sandpaper, the normal load signal fluctuates but its fluctuation amplitude is much smaller than its mean value when the gain is set at 0.001.This is shown by the black curve in Supplementary Fig. 6, confirming that the scanning is made under the condition that the normal load is kept approximately constant.
In the experiment, we choose the normal load N in the range of 200-600 nN so that the scanning probe is under a continuing stick-slip motion with irregular stick (and slip) periods.The scanning distance over the sandpaper is set at 100 µm for each scan so that one can collect an adequate amount of slip events to build up the probability density functions.The scanning speed U is kept in the low-speed regime so that individual slips can be fully resolved without any influence from other dynamic effects.With the chosen scanning speed U = 0.1 µm/s, each scan takes 16.7 minutes and collects 1 × 10 6 data points for the measured force trajectory F (x) (at the sampling rate 1 kHz).This is about the slowest scanning speed experimentally accessible for a mechanically stable AFM scan without much influence from noises resulting from mechanical instabilities and drifts of the system.We also vary the scanning speed U in the range 0.1-1 µm/s to verify the reproducibility of the experimental findings.
With our AFM setup, the response time t 0 of the scanning probe is t 0 ≃ √ m/k 0 , where m is the mass and k 0 is the static spring constant of the scanning probe.With m ≃ 1 × 10 −8 g estimated based on the dimension of the hanging beam and k 0 ≃ 0.48 N/m as measured for the 2D probe, we have t 0 ≃ 5 × 10 −6 s.For a typical slip distance λ 0 = 0.42 µm (see Supplementary Fig. 2(b)), we find a typical slip velocity λ 0 /t 0 ≃ 8.4 cm/s.Because the scanning speeds U used in the experiment are much smaller than the estimated slip velocity, our friction measurements are therefore all conducted in the low speed regime.Supplementary Fig. 7 shows a comparison of two repeated frictional force measurements at two different scanning speeds.It is seen that the two force curves overlap with each other in most positions, indicating that the measured F (x) and the resulting stick-slip events are not sensitive to the change of scanning speeds used in the experiment.There are a few locations where the two force curves show some deviations.These deviations are mostly associated with small-sized slips and can be attributed to the stochastic slip dynamics.As will be shown in Sec.II.F below, the obtained statistical properties of stick-slip friction at meso-scale are robust and do not change with the scanning speeds in the low speed regime.

II. SUPPLEMENTARY DISCUSSION ON DATA ANALYSIS A. Identification of slip events
Supplementary Fig. 8 shows a typical frictional force loop obtained when the 2D probe is pulled back and forth for a cycle.Before the scan starts, the 2D probe is in static contact with the sandpaper surface under a normal load N = 600 nN.As the probe is pulled to advance (→) (or recede (←)), the resulting static frictional force F increases (or decreases) quickly and linearly with the sandpaper displacement x and gives rise to a straight line on the left (or right) side of the force loop.When the pulling force reaches a critical value F c , the solid interface depins and the 2D probe begins its stick-slip motion, as shown by the top (or bottom) horizontal fluctuating force curve.
To identify the beginning and end of each slip event, we take a differential of the measured advancing frictional force F (x) and obtain the force differential curve dF (x) = F (x i ) − F (x i+1 ), where i and i + 1 are two adjacent sampling points of the measured F (x), as shown in Supplementary Fig. 9(a).Because the sampling rate is 1 kHz, the time step between the two adjacent sampling points, ∆t = 1 ms, is held constant.The corresponding spatial step of the displacement is ∆x = x i −x i+1 ≃ U ∆t.For U = 100 nm/s, we have ∆x = 0.1 nm.Because U is not held strictly constant (having small fluctuations), ∆x = 0.1 nm is true only on the average sense.As shown in Supplementary Fig. 9(b), individual slip events can be identified by the sharp spikes in the measured dF (x) curve.On average, a single slip (spike) lasts for about 5 ms (5 sampling points).
To discriminate the noise in the measured dF (x) curve, we set a threshold value ϵ, so that all the spikes in the measured dF (x) with an amplitude dF (x) > ϵ are counted as slip events.The threshold value ϵ is deter-

mined as follows. (i)
We first calculate the RMS value of dF (x) and remove those data points whose amplitude dF (x) is larger than the RMS value.(ii) We then repeat this procedure to the remaining data points for three times so that the final RMS value converges to a constant value, which is the intrinsic noise of the system made up mainly by the environmental vibration and electronic noise of the AFM.(iii) The threshold value ϵ is set as twice the final RMS value.The obtained threshold value ϵ is usually in the range of 0.2-0.5 nN, which varies slightly with the cantilever used and the environmental noise.Once the slip events are identified, the force at the beginning of each slip event is recorded as the maximal (depinning) force F c and the force drop between the beginning and end of each slip event is recorded as the force release δf associated with that slip event.

B. Determination of interfacial adhesion
Supplementary Fig. 10 shows four typical contactforce-versus-distance curves F (z) measured when the scanning probes approach and then recede from different substrates.The measured F (z) all increases linearly with the cantilever displacement z and the slope gives the spring constant of the cantilever in the normal direction.The measured force loops in the advancing and receding Supplementary Fig. 10.Measurements of interfacial adhesion.Plots show the measured vertical (contact) force as a function of vertical traveling distance of the cantilever when the scanning probe approaches and then recedes from different substrates.The black and red curves are obtained, respectively, for the quasi-1D probe and 2D probe over the sandpaper surface.The blue and green curves are obtained for the 2D probe over the gold-coated and bare silicon wafer surfaces, respectively.All the measurements are made at a constant speed of 2 µm/s.For clarity, the red, blue and green curves are shifted to the right by 1.5 µm (red), 2.9 µm (blue) and 4.4 µm (green), respectively.
directions overlap quite well with little hysteresis.It is seen from the black and red curves that when the scanning probe detaches from the sandpaper surface, there is a pull-off force resulting from the adhesion between the probe and sandpaper.This adhesion force is about the same for the quasi-1D and 2D probes over the sandpaper surface.By repeating the force measurements at 10 different locations on the same substrate, we obtain the mean value and the standard deviation of the adhesion force to be 21± 6 nN (for 1D probe/sandpaper) and 23±7 nN (for 2D probe/sandpaper), respectively.This adhesion force is small when compared with the normal load N used in the experiment, which is in the range of 200-600 nN.On the other hand, the adhesion force between the 2D probe and the smooth gold-coated silicon wafer surface as well as the bare silicon wafer surface increases significantly with the adhesion force equal to 59±4 nN (for 2D probe/gold-coated silicon wafer) and 125±2 nN (for 2D probe/silicon wafer), respectively.We, therefore, conclude that the reduction of adhesion between the scanning probe and sandpaper results mainly from the rough and chemically inert surface of the sandpaper.

C. Determination of the friction coefficient µc
Supplementary Fig. 11 shows how the mean value of the maximal lateral force F c needed to trigger the lo- cal slips changes with the normal load N .The relatively large error bars result from the large fluctuations of the maximal frictional force measured in the stick-slip regime.It is seen that the obtained ⟨F c ⟩ is linearly proportional to N , indicating that the Amontons' law of friction [11] still holds for stick-slip friction.The solid lines in Supplementary Fig. 11 show the linear fits, ⟨F c ⟩ = µ c N , to the data points with µ c = 0.55 ± 0.03 for the quasi-1D probe and µ c = 0.26 ± 0.02 for the 2D probe.

PDFs of the normalized maximal force fc
The generalized extreme value (GEV) distribution is given in Eq. ( 1) of the main text and is repeated here for ease of reference [12]: where f c is a random variable.In the general case in which f c is not normalized, Eq. ( 5) has three independent parameters: a location parameter µ, a scale parameter β and a shape parameter ξ.The location parameter µ and scale parameter β are approximately the mean and standard deviation of f c .The shape parameter ξ determines the tail shape of the GEV distribution, with ξ = 0 for an exponential tail, ξ > 0 for a heavy tail and ξ < 0 for a short or no tail.When f c is normalized with a zero mean and unity variance, as stated in Eq. ( 1) of the main text, the three parameters in Eq. ( 5) are bounded by the two normalization constraints, , where Γ(• • • ) is the Gamma function.In this case, Eq. ( 5) has only one independent parameter.By choosing ξ as the independent parameter, we have When ξ = 0, Eq. ( 5) can be simplified by taking the limit ξ → 0, i.e., ) , (6) where β 0 = √ 6/π and µ 0 = − √ 6γ/π with γ ≃ 0.577 being the Euler's constant.Equation (6) represents a special case of the GEV distribution (with ξ = 0) and is called the Gumbel distribution, which has a long exponential tail.Supplementary Figs.12(a) and 12(b) show the measured PDFs P (f c ) of the normalized maximal force f c for the quasi-1D probe and 2D probe, respectively, under different normal loads N .Here we define f c = (F c −⟨F c ⟩)/σ Fc with ⟨F c ⟩ and σ Fc being, respectively, the mean and standard deviation of F c .It is seen that for the same probe, the measured PDFs P (f c ) at different N all collapse onto a master curve, which can be well described by the GEV distribution in Eq. ( 5).The black solid lines in Supplementary Fig. 12 show the fits of Eq. ( 5) to the data points, from which we find ξ = −0.13± 0.06 for the quasi-1D probe and ξ = −0.03± 0.05 for the 2D probe.The error bar of the fitting parameters quoted here represents the confidence interval that we obtain from the fitting.The obtained values of ξ for the two probes are both close to zero, suggesting that the obtained PDFs P (f c ) are close to the Gumbel distribution in Eq. ( 6).Supplementary Fig. 12 thus suggests that the GEV (or Gumbel) distribution is a universal description of the depinning force for local slips, which is invariant with the normal load N in the stick-slip regime studied.

PDFs of the slip length δxs
Supplementary Figs.13(a) and 13(b) show the measured PDFs P (δx s ) of the slip length δx s for the quasi-1D probe and 2D probe, respectively, under different normal loads N .The slip length δx s is the force release δf normalized by the static spring constant k 0 , i.e., δx s = δf /k 0 .It is seen from Supplementary Fig. 13(a) that the measured PDFs P (δx s ) for the quasi-1D probe at different values of N all collapse onto a master curve, which can be well described by the power-law distribution for all values of δx s smaller than a cutoff slip length δx * s ≃ 800 nm.The fitted value of the power-law exponent is τ = 1.12 ± 0.10 (black solid line) and the error bar quoted here represents the confidence interval that we obtain from the fitting.Supplementary Fig. 13(a) thus suggests that the observed power-law distribution of the slip length δx s for the quasi-1D probe is invariant with the normal load N in the stick-slip regime studied.
It is seen from Supplementary Fig. 13(b) that the measured PDFs P (δx s ) for the 2D probe reveal some weak dependence on N for the two sets of data with N ≤ 400 nN.These two sets of data can still be fitted to Eq. ( 7) but the fitted value of τ decreases slightly with increasing N from τ = 0.92 at N = 250 nN to τ = 0.82 at N = 400 nN.Furthermore, the cutoff slip length δx * s increases gradually with N from δx * s ≃ 160 nm at N = 250 nN to δx * s ≃ 270 nm at N = 400 nN.For N ≥ 500 nN, the measured PDFs P (δx s ) reach to an asymptotic form and do not change with N any more.The two sets of data can be well described by a common power-law distribution, as shown in Eq. ( 7), with τ = 0.72 ± 0.10 for all values of δx s smaller than the cutoff slip length δx * s ≃ 400 nm.These results suggest that the observed power-law distribution of the slip length δx s for the 2D probe becomes invariant with the normal load N , when N is large enough so that the 2D probe is in an "optimal" contact with the sandpaper and can feel its entire roughness landscape.

PDFs of the local force gradient k ′
Supplementary Figs.14(a) and 14(b) show the measured PDFs P (k ′ ) of the local force gradient k ′ ≡ dF i (x s )/dx s for the quasi-1D probe and 2D probe, respectively, under different normal loads N .Here F i (x s ) is the random pinnning force field sampled by the solid interface and k ′ is extracted from the measured dynamic spring constant k using the equation k ′ = kk 0 /(k 0 − k) (see Eq. ( 13) below).In the plot, k ′ is normalized by the static spring constant k 0 , which is independent of N .It is seen that for the same probe, the measured PDFs P (k ′ ) at different values of N all collapse onto a master curve within the statistical error and can be described by the exponential distribution where the decay rate b is a fitting parameter.The red solid lines show the fits of Eq. ( 8) to the data points, from which we find b = 0.14 ± 0.02 (or ⟨k ′ ⟩/k 0 = 7.1) for both the quasi-1D probe and the 2D probe.

E. PDFs of fc, δxs and k at the high-load limit
To illustrate the different behaviors of stick-slip friction between the intermediate and high-load regimes, we consider the evolution of asperity configurations at the contact interface with increasing normal load N .Supplementary Fig. 15 shows a sketch of the contact geome-Supplementary Fig. 15.Evolution of contact geometry with increasing normal load N. The blue shaded area represents a deformable flat surface and the grey shaded area represents a rigid rough surface.
try between a deformable flat surface (e.g., the scanning probe with a Young's modulus of ∼ 3 MPa) and a rigid rough surface (e.g., the sandpaper with a Young's modulus of ∼ 170 GPa).In the low-load regime, the asperities (or microcontacts) at the interface are dilute and their number increases with the normal load N [13].When the asperities become dense enough, further increase of N introduces interactions and coalesce between nearby asperities.As a result, the number of asperities at the interface starts to decrease with N and consequently the incident rate of slips is reduced, as shown in Figs.1(e) and 1(f) of the main text.In fact, the effect illustrated in Supplementary Fig. 15 should apply to any normal load, so long as the rough surfaces has a broad range of asperity sizes and separations (see Sec. I.A).Under small normal loads, only those asperities with a very close separation can interact with each other.As the normal load increases, asperities with a relatively large separation can get in touch with each other due to their elastic and (or) plastic deformations.In other words, the coalescence effect is enhanced with an increasing N .
Assuming the actual contact area A r is ∼ 1% of the apparent contact area A a = 12 × 12 µm 2 = 1.44 × 10 −10 m 2 , we find the effective stress σ N exerted on the microcontacts at the normal load N = 2400 nN to be, σ N = N/A r ≃ 1.7 MPa.This value of σ N is close to the yield stress of the UV cured glue (estimated as the Young's modulus for a 100% strain).Supplementary Fig. 15 thus provides an argument, which extends the Greenwood and Williamson model [13] for interacting asperities in the high-load regime.The evolution of the contact geometry with an increasing N also implies a continuous change of the random pinning force field.Microscopically, because the local pinning sites in the high-load regime become stronger, multi-asperity slips where the entire scanning probe crosses over a number of asperities in a single slip event are more likely to occur, which also reduces the in- cident rate of slips.A similar effect of multi-asperity slips has been observed in nanotribology [14].As a result, both the effects of asperity coalescence and multi-asperity slips are possible reasons for the observed stick-slip behavior in the high-load regime.
The vertical deflection ∆z of the cantilever at the normal load N = 2400 nN can be estimated as, ∆z = N/k n ≃ 3.2 µm, where k n = 0.74 N/m is the spring constant of the cantilever in the normal direction.The cantilever used in the experiment has a rectangular shape and its length L = 300 µm.The resulting strain of the cantilever is ∆z/L ≃ 1.1%.For this small value of strain, vertical deformations of the cantilever are in the linear response regime and their nonlinear effect is negligible.The measurement of the frictional force, on the other hand, uses the lateral deflection ∆x of the cantilever, which is independent of ∆z.It is seen from Fig. 1(f) in the main text that the measured frictional force F (x) in the stick state changes linearly with the traveling distance x, indicating that the lateral deformations of the scanning probe are also in the linear response regime.Slight deviations of the measured F (x) from its linear dependence on x are observed when the measured F (x) reaches its maximal value before slips take place.We believe that this is caused by the partial slip (or plastic deformation) of the local contacts but not by the nonlinear effect of the hanging-beam probe.This is because a similar behavior is also observed when the value of N is in the intermediate range (200-600 nN), in which the amplitude of the measured F (x) is reduced by more than 4 times, and also in other experiments as reported in Ref. [15].Supplementary Figs.16(a), 16(b) and 16(c) show, respectively, the PDFs of the three quantities, f c , δx s and k ′ , measured at the high load limit, N = 2400 nN.The corresponding force trajectories are shown by the red curve in Fig. 1(e) of the main text.The measured P (f c ) in Supplementary Fig. 16(a) has a nearly symmetric shape, which can be described equally well by either the GEV distribution in Eq. 5 with a large value of ξ ≃ −0.26 (red solid line) or a (Gaussian) normal distribution (black solid line).The measured P (f c ) at the high load limit is very different from that obtained at the intermediate values of N , which follows the GEV distribution with a small value of ξ and a long tail for the positive values of f c .Note that the obtained values of f c in Supplementary Fig. 16(a) are distributed in a narrower range [-2, +2] compared with that measured at the intermediate values of N , suggesting that the variation range of the resulting f c at the high load limit is decreased.
The measured P (δx s ) in Supplementary Fig. 16(b) does not show a prominent power-law distribution, as that obtained at the intermediate values of N .It reveals a noisy decay for small values of δx s (< 0.1 µm), followed by a plateau region when δx s is in the range between 0.1 µm and 1 µm.This behavior of P (δx s ) suggests that at the high load limit, slips tend to have larger sizes compared with those obtained at the intermediate values of N .
The measured P (k ′ ) in Supplementary Fig. 16(c) does not follow the exponential distribution, as that obtained at the intermediate values of N .It reveals a plateau region for small values of k ′ /k 0 (< 8), followed by a rapid  The top panels (a,b,c) are measured using the quasi-1D probe at a normal load N = 500 nN and the bottom panels (d,e,f) are obtained using the 2D probe at a normal load N = 500 nN.The solid lines in (a) and (d) are the fits of the GEV distribution in Eq. ( 5) to the data with (a) ξ = −0.13± 0.06 and (d) ξ = −0.03± 0.05.The solid lines in (b) and (e) are the fits of the power-law distribution in Eq. ( 7) to the data with (b) τ = 1.12 ± 0.10 and (e) τ = 0.72 ± 0.10.The solid lines in (c) and (f) are the fits of the exponential distribution in Eq. ( 8) to the data with (c) b = 0.14 ± 0.02 and (f) b = 0.14 ± 0.02.
decay for larger values of k ′ /k 0 (> 8).This behavior of P (k ′ ) suggests that the scanning probe experiences a stronger pinning field at the high load limit and thus produces a heavier distribution than the exponential one.This is consistent with the observation that the measured PDF P (k) (inset of Supplementary Fig. 16(c)) of the dynamic spring constant k has a higher weight at large k (≥ 0.6k 0 ) compared with that measured at the intermediate values of N .
The above results indicate that the stick-slip friction at large values of N does not follow the three statistical laws as observed at the intermediate values of N .This is caused by a large reduction both in numbers and in strengths of the asperities at the contact when it is under a large normal load.This finding further suggests that the stick-slip friction under the intermediate normal load is at a "critical state", in which the contact between the two solid surfaces has an adequate number of random asperities and each individual asperity can locally generate a unique pinning to the moving interface.Under this condition, the scanning probe explores the entire landscape of the pinning force field at the contact with the rough sandpaper surface, so that the random stick-slip friction is characterized by the three statistical laws as reported in the main text.

F. Scanning speed effects on the frictional force statistics
Supplementary Fig. 17 shows the measured PDFs of the normalized maximal force f c (a,d), slip length δx s (b,e) and local force gradient k ′ (c,f) for different scanning speeds U .It is seen that the measured PDFs at different values of U all collapse onto a single curve, suggesting that the functional form of the measured P (f c ), P (δx s ) and P (k ′ ) is invariant with the scanning speed U studied.This conclusion applies to all the friction measurements conducted by using both the quasi-1D and 2D probes.The measured P (f c ), P (δx s ) and P (k ′ ) are well described, respectively, by the GEV distribution in Eq. ( 5), the power-law distribution in Eq. ( 7) and the exponential distribution in Eq. ( 8) with the same fitting parameters as those reported in the main text.The results shown in Supplementary Fig. 17 thus demonstrate that the observed statistical laws of stick-slip friction at meso-scale are robust and do not change with the scanning speed U in the low speed regime.To further understand the low-load behavior of the frictional force at mesoscale, we conduct additional measurements of the frictional force in the small load regime (N < N m with N m = 124 nN, see the main text).Supplementary Fig. 18(a) shows the frictional force loops under three different normal loads in the small load regime.At N = 100 nN, the force trajectory shows continuous fluctuations with a few abrupt changes (slips) occasionally.This kind of force trajectories are considered as smooth sliding, because the scanning probe only feels a weak pinning field with a finite number of predominantly small asperities, which are too small to cause stick-slip instabilities.As N decreases, the number of small asperities decreases (see Supplementary Fig. 15) and so does the magnitude of the measured force fluctuations.The residual rare slip events become fewer and eventually disappear at N = 10 nN.Note that the hysteresis between the mean frictional forces in the advance direction and in the receding direction also decreases with decreasing N .Such a reduction of friction-induced hysteresis with decreasing normal load was also observed in atomic friction at the single contact limit [16].Supplementary Fig. 18(b) shows how the measured frictional force loops change with the scanning speed U when the normal load is held at a small constant value of N = 20 nN.It is seen that the measured frictional force loops do not reveal any significant speed-dependence in the small load regime.This result is in agreement with a general trend that was observed for rough surfaces [15,17] that the dynamic friction coefficient in steady sliding is roughly a constant or decreases slightly with increasing U when U is small, with a nonzero quasi-static limit as U → 0. This general trend is also confirmed in the high-load regime, as shown in Supplementary Fig. 17.The result shown in Supplementary Fig. 18(b), however, does not agree with those obtained in atomic friction at the single contact limit, in which the frictional force goes to zero as U → 0 [16,18].
For weak pinning fields, such as those in atomic friction with a characteristic energy barrier E b and lattice constant λ, the stick-slip instability gives rise to a nonzero frictional force of the order of E b /λ at the low sliding speed limit [18,19].For even weaker pinning fields (or smaller normal load N ) where the stick-slip instability does not occur, the sliding tip can follow the corrugated surface landscape and the effect of the energy barrier E b may be averaged out with a near-zero frictional force at the limit of low sliding speed and single elastic contact [16,18].The situation is changed when the apparent contact area contains a finite number of micro-contacts or asperities, which are randomly distributed both in space and in size.In this case, the effect of the energy barriers associated with different asperities cannot be averaged out coherently, and thus gives rise to a non-zero frictional force at the quasi-static limit as U → 0.
As indicated in Eq. ( 14) below, the lateral force measured by AFM is a sum of the viscous drag force and the pinning (frictional) force.The former is a linear function of U , whereas the latter usually does not change much with U , as discussed above.In most dry friction experiments, the viscous drag force at the low speed limit is so small compared with the frictional force and thus is neglected [20,21].The results shown in Supplementary Fig. 18(b) thus suggest that force fluctuations in the small load (or steady sliding) regime result primarily from the low-load friction at the mesoscale contact, which contains a finite number of asperities and gives rise to a weak pinning force that has a non-zero value even at the low speed limit.

III. SUPPLEMENTARY DISCUSSION ON THEORETICAL MODEL
A. Frictional force and effective spring constant in the stick state The total frictional force at a solid interface can be expressed in terms of the asperity-induced heterogeneous tion between x C and x D , say x E , and then undergoes a damped oscillation and gradually drifts toward x C , until the remaining kinetic energy is fully dissipated.Because the slip occurs rapidly, the pulling point remains at x (2) 0 during the slip and the displacement of the CoM (or slip length) is thus given by δx s = x E − x B .Because of friction at the solid interface, the actual stick-slip motion is always influenced by some level of damping.As will be shown below, the use of the relation, δx s = δf /k 0 , to estimate the slip length, which is valid for over-damped systems, is indeed a good approximation for our system with an accuracy better than 90%.
By taking a differential of the stick-state condition, F e (x 0 ; x s ) = F i (x s ), with respect to x s , we have dx s = dx 0 /(1 + k ′ /k 0 ), where k ′ = dF i (x s )/dx s is the local (upward) force gradient (or spring constant) of the pining force field F i (x s ), as indicated in Supplementary Fig. 19.This relation suggests that the CoM is pinned only partially in the stick state and its movement dx s is a factor of 1/(1 + k ′ /k 0 ) smaller than the displacement δx 0 of the pulling point.Furthermore, the instantaneous velocity dx s /dt of the CoM is smaller than U = dx 0 /dt when the CoM moves upwards (k ′ > 0) along F i (x s ) and is larger than U when it moves downward (k ′ < 0).This incoherent motion between the CoM and the pulling spring gives rise to a fluctuating dynamic spring constant k defined as k ≡ dF i (x s )/dx 0 .By taking a derivative of F e (x 0 ; x s ) = F i (x s ) with respect to x 0 , we have Equation (13) indicates that the measured k is an effective spring constant of two springs connected in series, wherein one has a spring constant k 0 and the other has k ′ .The onset condition for slip (i.e., k = k 0 ) can be experimentally verified by examining the steady-state condition that the total displacement ∆x s of the CoM (of the scanning probe) over a long run is equal to the total displacement ∆x 0 of the pulling point.The latter is just the scanning distance of the AFM scanner, i.e., ∆x 0 = 500 µm, which can be directly read out from the AFM controller.The former contains two parts, i.e., ∆x s = ∆x s is the total displacement during the partial slip in the stick state.With these definitions, we have ∆x , where (δx s ) i and δf i are, respectively, the slip displacement and force drop associated with the ith slip event, and ∆x , where (δx 0 ) i and k ′ i are, respectively, the displacement of the pulling point and the local force gradient of the pinning force during the ith stick event.From the force trajectories obtained at U = 100 nm/s and under different normal loads N , we calculate the values of ∆x s , and the final results are summarized in Supplenmentary Table 1.It is seen that the total displacements ∆x s of the CoM obtained under different experimental conditions are all very close to the scanning distance ∆x 0 (with over 91% agreement).This result thus confirms that the relation δx s ≃ δf /k 0 , which is used to calculate ∆x s , indeed gives a good estimate of the slip length from the measured force drop.The small difference between ∆x s and ∆x 0 (< 9%) may be attributed to a small underestimation made by using the relation δx s ≃ δf /k 0 , as discussed above, and/or the nonlinear effect of very large slip events.

B. An under-damped spring-block model for the slip dynamics
We now propose an under-damped spring-block model, as illustrated in Supplementary Fig. 19, to explain the observed power-law distribution of the slip length δx s .The equation of motion of the CoM position x(t) can be written as (here we use x instead of x s to simplify the notation) The first term on the left hand side of Eq. ( 14) is the inertial term with m being the effective mass of the scanning probe.The second term is a damping term with γ being the damping coefficient.This term is used to describe rapid dissipations in solid friction.This term was first adopted in the original Prandtl and Tomlinson (PT) model [16,22] and later in many of the follow-up models [15,23,24].On the right hand side of Eq. ( 14), the first term is the elastic pulling force with k 0 being the static spring constant.The pulling force is exerted by the cantilever and hanging beam moving at a constant speed U .The second term F i (x) (> 0) is the asperity-induced random pinning (or frictional) force field acting on the CoM.Note that F i (x) is intrinsically dissipative and its mean value F 0 = ⟨F i (x)⟩ > 0 is responsible for keeping the system at a steady state.The constant frictional force F 0 only gives rise to an additional steady-state stretching F 0 /k 0 to the pulling spring without changing the fluctuation dynamics.The fluctuating force F i (x) − F 0 is responsible for the dynamics of the scanning probe.It is assumed that where D is a measure of the fluctuation amplitude of F i (x).The Brownian correlated random force field has been used previously in the Alessandro-Beatrice-Bertotti-Montorsi (ABBM) model for the over-damped avalanche dynamics [25,26].
As mentioned in the main text, Eq. ( 14) is a further extension of the PT model, which has been widely used in the study of dry friction and nanotribology [16,22].The PT model was considered as the most successful and influential "minimalistic" model in rationalizing the wealth of nanoscale and mesoscale friction data produced over  s , which varies between the quasi-1D probe and 2D probe.The scanning distance of the AFM scanner is ∆x0 = 500 µm, which is the same for both the quasi-1D and 2D probes.
the last decades [18,22,27].A number of experiments have been conducted and demonstrated the characteristic underdamped features of dry friction.For example, attenuated oscillations after a slip were observed in the friction-force traces obtained using a surface force apparatus with two mica surfaces sliding across a solidified molecular film [28].Overshoots of a single slip over multiple slip-lengths were found in the lateral force trajectories obtained when a sharp AFM tip slides over a crystalline solid surface [14].These measurements clearly demonstrated that dry friction in the single-asperity limit is underdamped and thus the PT model is deemed appropriate to describe many fundamental properties of dry friction.
The present experiment broadens the study of dry friction into the regime of multi-asperity dynamics and so does our model, which provides an intermediate level description bridging the gap between the microscopic behavior of individual slips and the macroscopic laws of solid friction.Because the decay rate that determines how long the inertia effect will last scales as γ/m, which is independent of the system size (or mass), dry friction in the multi-asperity regime will remain underdamped as that in the single-asperity regime.As a result, the predictions of the ABBM model and its variations, which are valid only for the overdamped systems, will not be applicable to dry friction at mesoscale.
In recent years, many attempts have been made to include some new features of the random force field into the PT model, such as adding a Gaussian white noise to the PT model to account for the thermal activation effect [15], assuming F i (x) results from a fractal potential [29], or assuming F i (x) to follow the Amontons' law of friction, F i (x) = µN , with the frictional coefficient µ following the (empirical) rate (or speed-dependent) and state (or time-dependent) friction law [23,24].These models studied the transition boundary between the steady sliding and stick-slip motion, as predicted by the PT model, and provided a more realistic but complex framework for the pinning force field.However, to our knowledge, these modified PT models have not reported any result on the slip length distribution.
When a sharp AFM tip scans over a fractal-like rough surface, as proposed in Ref. [29], the tip can feel the entire roughness landscape of the underlying surface.In this case, the resulting slip lengths may follow a powerlaw distribution, which is linked to the statistics of the surface roughness or the pinning force field.Because it is solely determined by the self-similar nature of the surface roughness, the power-law distribution in this singleasperity regime does not have any connection with the avalanche dynamics.In this work, we demonstrate that the stick-slip friction in mesoscale is an emergent (or selforganized) phenomenon, which is realized in the multipleasperity regime.In this case, the power-law distribution of slip length δx s arises from the spontaneous response of multiple asperities at the contact to the stick-slip instabilities and is not necessarily linked to the self-similar nature of the surface roughness.Instead of assuming the total pinning force F i (x) in Eq. ( 14) has a specific fractal form, our model adopts a more general assumption that the total force F i (x) resulting from the multiple-asperity pinning is Brownian correlated.
In fact, the assumption that F i (x) is Brownian correlated is a general feature of the random pinning force field.To show this, we rewrite Eq. ( 9) as [σ(x, y; x s )dy]dx = ∫ Lx h(x; x s )dx, (15) where h(x; x s ) = ∫ Ly σ(x, y; x s )dy is an effective interfacial tension.Because both σ(x, y; x s ) and F i (x s ) are real physical quantities, they must have a finite variance.Consequently, the integrand h(x; x s ) in Eq. ( 15) will also have a finite variance σ 2 h , which is related to the frictional force variance σ 2 F by σ 2 F = σ 2 h L x .In fact, this statement is true even for fractal-like rough surfaces, because the self-similar feature of the actual solid surfaces always has two cutoff lengths.The lower cutoff is determined by the smallest asperity size (e.g., the grain size of the sandpaper), whereas the higher cutoff must be bounded by the system size.
The last equality in Eq. ( 15) states that the total pinning force F i (x s ) is a running sum of h(x; x s ) along the moving direction.According to the running average (or sum) theorem [30], the auto-correlation function c(λ) = ⟨F i (x s )F i (x s + λ)⟩ of F i (x s ) has the form where L x − λ is the overlap region between F i (x s ) and F i (x s + λ) (for small values of λ/L x ).The mean-squared force-difference ⟨∆F 2 i (λ)⟩ = ⟨|F i (x s + λ) − F i (x s )| 2 ⟩ then takes the form where D = σ 2 F /L x = σ 2 h .Differentiating both sides of Eq. ( 14) with respect to x and using d/dx = 1 v d dt , where v (≡ dx/dt) is the velocity of the CoM, we have where w(x) = dF i (x)/dx has a zero mean and obeys the relation ⟨w(x)w(x ′ )⟩ = 2Dδ(x−x ′ ).Multiplying Eq. ( 18) by v, we obtain the equation of the slip velocity v(t), where η(t) = √ vw is a Gaussian white noise, having a zero mean and a correlation ⟨η(t)η(t ′ )⟩ = 2Dδ(t − t ′ ).Now making a variable replacement with v = V U and t = T √ m/k 0 , we obtain The new noise ξ(T ) has a correlation ⟨ξ(T )ξ(T ′ )⟩ = (2D √ m/k 0 )δ(T − T ′ ), which can be obtained by using the correlation function ⟨η

C. Numerical results of the under-damped spring-block model
Equation ( 21) has two parameters, the normalized damping coefficient γ ′ and the normalized "diffusion coefficient" D ′ of the random pinning force field.In the absence of the random force field, i.e., when D ′ = 0, V oscillates and converges to 1 for γ ′ ≤ 2, indicating that the system is in the under-damped regime.For γ ′ > 2, V monotonically decays to 1, suggesting the system is over-damped.On the other hand, when D ′ ≫ 1, Eq. (21) suggests that the motion of the CoM is dominated by the pinning-depinning dynamics in both the under-damped and over-damped regimes.
Equation ( 21) can be solved numerically using the Euler-Maruyama method, following the protocol adopted in Ref. [26].A slip starts at a zero acceleration ( dV dT (T = 0) = 0) and a very small positive initial velocity V 0 ≪ 1 and stops when V first becomes negative after a period T s .Since the velocity should always be positive, a new slip starts on the same initial condition but with a different realization of the noise.The total displacement during a slip is defined as the slip length, δx s ≡ ∫ Ts 0 V dT .The initial velocity V 0 used here only affects the lower cut-off of the power-law regime of the slip-size distribution P (δx s ).It is found that a smaller value of V 0 gives rise to a smaller value of the lower cut-off.We chose V 0 = 10 −4 in this work in order to have a large powerlaw range.Supplementary Fig. 20 shows four examples of the velocity time series V (T ) during a slip, which are obtained at γ ′ = 2 and D ′ = 500.It is seen that these slip events have different slip durations T s and different maximum velocities.The velocity trajectories show an overall acceleration phase followed by a deceleration phase, but the actual shape of V (T ) varies for different slips, which is indeed a manifestation of the stochastic nature of avalanches.By carefully examining the individual velocity trajectories, we find that neither the velocity V (T ) nor the slip duration T s remain constant during the slips.Supplementary Fig. 21 show the PDFs P (δx s ) of the numerically calculated slip length δx s for different values of D ′ , when the system is in the under-damped regime with γ ′ = 1.It is seen that in the strong pinning limit D ′ ≫ 1 (e.g., D ′ = 500), the obtained P (δx s ) follows a power-law distribution, P (δx s ) ∼ δx −ϵ s , over a wide range of δx s for seven decades.Note that the power-law distribution of slip length δx s arises from the dynamics of stick-slip instabilities and is not necessarily linked to a specific feature of the surface roughness.There is only a general requirement for the rough surface under study that the corresponding pinning force field F i (x) is random enough and Brownian-correlated.Therefore, the stick-slip friction in mesoscale is an emergent (or selforganized) phenomenon, which is realized in the multipleasperity regime.
As D ′ decreases, the power-law range is reduced mainly Analysis 6 A. Identification of slip events 6 B. Determination of interfacial adhesion 7 C. Determination of the friction coefficient µ c 7 D. Normal load effects on the frictional force statistics 8 1. PDFs of the normalized maximal force f c 8 2. PDFs of the slip length δx s 9 3. PDFs of the local force gradient k ′ 9 E. PDFs of f c , δx s and k at the high-load limit 10 F. Scanning speed effects on the frictional force statistics 12 G.Frictional force fluctuations at mesoscale in the small load regime (N < N m ) 13III.Supplementary Discussion on Theoretical Model 13 A. Frictional force and effective spring constant in the stick state 13 B.An under-damped spring-block model for the slip dynamics 15 C. Numerical results of the under-damped spring-block model 17 D. Steady-state velocity time series of the underdamped spring block model 18 E. Numerical results of the over-damped spring-block model 19 Supplementary References 19

1 .
Surface morphology of the sandpaper.(a) An AFM topographical image of the ultrafine sandpaper surface over an area of 50 µm × 50 µm.The horizontal scale bar is 10 µm and the vertical grey scale bar shows the surface height of the sandpaper.(b) Variations of the surface roughness height h(x) along the red line marked in (a).

3 .
Surface morphology of the gold-coated silicon wafer.The AFM topographical image has an dimension of 10 µm×10 µm, where the vertical grey scale indicates the surface height.The RMS roughness of the gold-coated wafer surface is 1.0 nm over the whole scanning area.

Supplementary Fig. 4 .
Calibration of lateral force sensitivity of the scanning probe.(a) Calibration of a glass fiber, which is deformed under a normal force F applied by an AFM cantilever of known spring constant.(b) Calibration of the scanning probe using the glass fiber as a force sensor.The glass fiber is deformed under a lateral force F applied by the hanging beam probe.The red arrows indicate the direction of the applied forces.
MPa (with E ≃ 3 MPa and the Poisson's ratio ν = 0.5 for incompressible isotropic materials).Also with A max r = A a ≃ 100 µm 2 and ℓ ≃ 90 µm, we have k max t ≃ 3.3 N/m.With the estimated values of k c , k b , and k max t , we obtain the maximum effective lateral spring constant k max e ≃ 1.3 N/m using Eq. (

Supplementary Fig. 6 .
Influence of deflection feedback gain on the frictional force measurements.The two measurements are made at feedback gain values of 0.001 (red curve) and 0.1 (blue curve) under the same normal load N = 400 nN (mean value) and same scanning speed U = 100 nm/s.For clarity, the origin of the blue curve is shifted downwards by 250 nN.The circled portions of the blue curve show the feedback-induced instabilities.The black curve shows fluctuations of the normal force curve recorded simultaneously with the red curve.

Supplementary Fig. 7 .
Frictional force measurements at different scanning speeds.The two measurements are made at a scanning speed U = 0.1 µm/s (black curve) and U = 1 µm/s (red curve) when the 2D scanning probe slides over the same region on the sandpaper and under the same normal load N = 600 nN.

8 .
An example of the measured frictional force loop.The measurement is made when the 2D probe is pulled back and forth for a cycle under the normal load N = 600 nN and at scanning speed U = 100 nm/s.

Supplementary Fig. 9 .
Identification of slip events from a measured force trajectory.(a) A zoom-in plot of the measured advancing frictional force F (x) in Supplementary Fig. 8.(b) Corresponding force differential dF (x) between two adjacent sampling points of the measured F (x) in (a).

Supplementary Fig. 11 .
Measurements of the friction coefficient µc.The two sets of data are measured for the quasi-1D probe (black circles) and 2D probe (red triangles) at a scanning speed U = 100 nm/s.The error bars show the standard deviation of the measured Fc.The solid lines show the linear fits, ⟨Fc⟩ = µcN , to the data points with µc = 0.55 ± 0.03 for the quasi-1D probe and µc = 0.26 ± 0.02 for the 2D probe.

Supplementary Fig. 12 .
PDFs of the normalized maximal force under different normal loads.The measurements are made for the quasi-1D probe (a) and 2D probe (b) at the same scanning speed U = 100 nm/s.The black solid lines show the fits of the GEV distribution (Eq.(5)) to all the data points with (a) ξ = −0.13±0.06 and (b) ξ = −0.03±0.05.

Supplementary Fig. 13 .
PDFs of the slip length under different normal loads.The measurements are made for the quasi-1D probe (a) and 2D probe (b) at the same scanning speed U = 100 nm/s.The black solid lines show the powerlaw fits, P (δxs) ∼ (δxs) −τ , to the data points with τ = 1.12± 0.10 (a) and τ = 0.72 ± 0.10 (b).The power-law fit in (a) is made to all the four sets of data points.The power-law fit in (b) is made only to the two sets of data points with N = 500 nN and N = 600 nN.

Supplementary Fig. 14 .
PDFs of the local force gradient under different normal loads.The value of k ′ is extracted from the measured dynamic spring constant k using Eq.(13) and is normalized by the static spring constant k0.The measurements are made for the quasi-1D probe (a) and 2D probe (b) at the same scanning speed U = 100 nm/s.The red solid lines show the fits of Eq. (8) to the data points with (a) b = 0.14 ± 0.02 and (b) b = 0.14 ± 0.02.

Supplementary Fig. 16 .
Statistical properties of the measured stick-slip friction at the high load limit.(a) Measured PDF P (fc) of the normalized maximal force fc.The red solid line shows a fit of the GEV distribution (Eq.(5)) to the data points with ξ ≃ −0.26.The data can also be described by a (Gaussian) normal distribution (black solid line).(b) Measured PDF P (δxs) of the slip length δxs.(c) Measured PDF P (k ′ ) of the local force gradient k ′ .The value of k ′ is extracted from the measured dynamic spring constant k using Eq.(13) and is normalized by the static spring constant k0.Inset shows the measured PDF P (k) of the normalized dynamic spring constant k/k0.All the measurements are made using the 2D probe at the normal load N = 2400 nN and scanning speed U = 100 nm/s.

20 .
Examples of the velocity trajectory during a slip.The four velocity time series V (T ) are obtained from the numerical simulation of Eq. (21) with γ ′ = 2 and D ′ = 500.

Table 1 .
Obtained values of the total displacement ∆x during the partial slip in the stick state, which is independent of N .The total displacement of the CoM is ∆xs = ∆x